Results

Transferring Data from Slicer to R

# Load necessary libraries #
library(devtools)
library(SlicerMorphR)
library(rgl)
library(geomorph)
library(readxl)
library(ks)

#  Set your working directory
setwd("INSERT_PATH_HERE")

##############################################
#   Load SlicerMorph Data and Preprocessing  #
##############################################

# Choose the Log File from Slicer
# (Ensure a line break is added at the end of the .log file before loading)

# Choose the SlicerMorph .log file
SM.log.file <- file.choose()  
SM.log <- parser(SM.log.file, forceLPS = TRUE)

# Extract Landmark Data
Data <- SM.log$LM

# Create geomorph dataframe
datagdf <- geomorph.data.frame(landmarks = Data)  

# Perform Procrustes Superimposition
Pcoords <- gpagen(datagdf$landmarks)

# Save Procrustes Coordinates
save(Pcoords, file = "filename.bin")

##############################################
#          Semi-Landmark Sliding            #
##############################################

# Extract Semi-Landmarks Data
Data2 <- SM.log$semiLMs

# Sliding using **bending energy minimization**
PcoordsSlid <- gpagen(A = Data, surfaces = as.numeric(Data2), ProcD = FALSE, print.progress = TRUE)
save(PcoordsSlid, file = "filenameSlidLandmarks.bin")


# Extract and Save Centroid Size (Csize)
write.table(PcoordsSlid$Csize, file = "FileName1.tsv", sep = "\t")


##############################################
#       Merge with Master Dataset            #
##############################################

# Load Pcoords Lists
PcoordsList <- read.table(file = 'FileName1.tsv', sep = '\t', header = TRUE)

# Load Master Dataset
Masterfile <- read_excel("path/to/YourMasterfile.xlsx")

# Merge Pcoords with Master Dataset
Table_SSL <- merge(PcoordsList, Masterfile, by = "LabID", sort = FALSE)

# Save Merged Tables
write.csv(Table_SSL, file = "FileName2.csv", row.names = FALSE)


##############################################
#   Create and Save Final Analysis Dataset   #
##############################################

# SSL Dataset

spec <- read.csv("FileName2.csv", header = TRUE)
YourName.dt <- list(gpa.sh = PcoordsSlid$coords, cs = PcoordsSlid$Csize, spec = spec)
save(YourName.dt, file = "filename.dt.bin")


##############################################
#            Symmetry Analysis               #
##############################################

## Define Landmark Pairs for Symmetry Analysis - Replace by your own paired landmarks
lm.pairs <- matrix(c(1, 2, 4, 5, 8, 9, 11:20, 22:85, 87:128, 130:177, 179, 180, 182, 183, 185, 186, 188:301, 304:399, 401:412, 414:427, 429:446, 448:461, 465:600, 602:625, 627:652, 656:665, 667:694, 696:701), ncol = 2, byrow = T)


## Perform Symmetry Analysis
asym <- bilat.symmetry(PcoordsSlid$coords, ind = dimnames(PcoordsSlid$coords)[[3]], 
                           object.sym = TRUE, land.pairs = lm.pairs, iter = 9)
symm.sh <- asym$symm.shape


# Save Symmetry Data

YourNewName.dt <- list(gpa.sh = PcoordsSlid$coords, cs = PcoordsSlid$Csize, 
                                   symm.sh = symm.sh, spec = spec)
save(YourNewName.dt, file = "YourNewfilename.dt.bin")

PCA

# Load necessary libraries
library(geomorph)
library(readr)
library(ks)
library(rgl)

# Set working directory (modify as needed)
setwd("INSERT_PATH_HERE")

# Load dataset
load("YourNewfilename.dt.bin")

######################################################################
#                          PCA Analysis                              #
######################################################################

# Perform Principal Component Analysis (PCA)
YourNamePCA <- gm.prcomp(YourNewfileName.dt$symm.sh)

# Access eigenvalues
eigenvalues <- YourNamePCA$sdev^2

# Calculate proportion of variance explained by each principal component
variance_explained <- eigenvalues / sum(eigenvalues)

# Print proportion of variance explained by each principal component
print(variance_explained)

######################################################################
#                          3D PCA Plot                               #
######################################################################

# Plot - Replace x by your number of groups in the categorical variable
cols <- rainbow(x)
plot3d(YourNamePCA$x[,1:3], col = cols[as.factor(YourNewfileName.dt$spec$YourGroup)], size = 20) 

# Add legend 
legend3d("topright", legend = unique(as.factor(YourNewfileName.dt$spec$YourGroup)), col=c(unique(cols[as.factor(YourNewfileName.dt$spec$YourGroup)])), pch = 10)

# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8) 


######################################################################
#                        3D Cloud Plot                               #
######################################################################

# Extract PCA scores
Score <- YourNamePCA$x[,1:3]
Labels <- factor(x=YourNamePCA$YourGroup) 

# Kernel density estimation for group classification
Hscv1.Bd <- Hkda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), bw = "scv", pre = "sphere")

# Perform kernel discriminant analysis
Tt.kda3.Bd <- kda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), Hs = Hscv1.Bd)

# Plot 3D Cloud
plot(Tt.kda3.Bd, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")

# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8) 

# Perform PERMANOVA IN PAST
# Export the Pc Scores #
write.table(YourNamePCA$x, file="FileName3.tsv", sep="\t") # tsv because of space


######################################################################
#                          Plot 3D Coastal Only                      #
######################################################################

# Filter out group of specimens (e.g., filter out offshore)

# Identify indices for specimens we want to keep
Indice <- YourNewfileName.dt$spec$YourGroup != "offshore"

# Subset the dataset based on offshore specimens
YourNewfileName.dt$symm.sh <- YourNewfileName.dt$symm.sh[, , Indice, drop = FALSE]
YourNewfileName.dt$cs <- YourNewfileName.dt$cs[Indice, drop = FALSE]
YourNewfileName.dt$spec <- YourNewfileName.dt$spec[Indice, ]

# Create Geomorph Data Frame (GDF) without offshore specimens
gdfNoOffshore <- geomorph.data.frame(
  shape = YourNewfileName.dt$symm.sh,
  ecotype = YourNewfileName.dt$spec$YourGroup,
  ind = YourNewfileName.dt$spec$LabID,
  Csize = YourNewfileName.dt$cs
)

# Perform PCA
YourNamePCA <- gm.prcomp(YourNewfileName.dt$symm.sh)

######################################################################
#                          3D PCA Plot                               #
######################################################################

# Replace x by your number of group in the categorical variable
cols <- rainbow(x)
plot3d(YourNamePCA$x[,1:3], col = cols[as.factor(YourNewfileName.dt$spec$YourGroup)], size = 20) 

# Add a legend
legend3d("topright", legend = unique(as.factor(YourNewfileName.dt$spec$YourGroup)), col=c(unique(cols[as.factor(YourNewfileName.dt$spec$YourGroup)])), pch = 10)

# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8)

######################################################################
#                        3D Cloud Plot                               #
######################################################################

# Extract PCA scores
Score <- YourNamePCA$x[,1:3]
Labels <- factor(x=YourNamePCA$YourGroup) 

# Kernel density estimation for group classification
Hscv1.Bd <- Hkda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), bw = "scv", pre = "sphere")

# Perform kernel discriminant analysis
Tt.kda3.Bd <- kda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), Hs = Hscv1.Bd)

# Plot 3D plot
plot(Tt.kda3.Bd, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")

# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8)


######################################################################
#              Plot lollipop graphs for PCA components               #
######################################################################

PC1 <- plotRefToTarget(
  YourNamePCA$shapes$shapes.comp1$min, 
 YourNamePCA$shapes$shapes.comp1$max, 
  method = "vector", 
  mag = 1, 
  gridPars = gridPar(pt.size = 0.25, pt.bg = "red")
)

PC2 <- plotRefToTarget(
  YourNamePCA$shapes$shapes.comp2$min, 
 YourNamePCA$shapes$shapes.comp2$max, 
  method = "vector", 
  mag = 1, 
  gridPars = gridPar(pt.size = 0.25, pt.bg = "blue")
)

PC3 <- plotRefToTarget(
  YourNamePCA$shapes$shapes.comp3$min, 
  YourNamePCA$shapes$shapes.comp3$max, 
  method = "vector", 
  mag = 1, 
  gridPars = gridPar(pt.size = 0.25, pt.bg = "green")
)

######################################################################
#              Export PCA scores for external analysis               #
######################################################################

# Export PCA scores for external analysis
write.table(YourNamePCA$x, file = "FileName4.tsv", sep = "\t")
A local Image
A local Image

Allometry Analysis

# Set working directory
setwd("INSERT_PATH_HERE")

# Load necessary libraries
library(geomorph)

# Load data
load("filenameSlidLandmarks.bin")
load("YourNewfilename.dt.bin")


###########################################
# HYPOTHESIS 1: ECOLOGICAL ALLOMETRY TEST #
###########################################


## Create Geomorph Data Frame (GDF) ##
gdf <- geomorph.data.frame(
  shape = filenameSlidLandmarks$coords, 
  YourGroup = filename.dt$spec$YourGroup, 
  ind = filename.dt$spec$LabID, 
  Csize = filenameSlidLandmarks$Csize
)

# Reduced (Null Hypothesis)
YourName_Reduced <- procD.lm(shape ~ log(Csize) + YourGroup, data = gdf, SS.type = "I", iter = 999)
anova(YourName_Reduced)

# Full Model (Alternative Hypothesis)
YourName_Full <- procD.lm(shape ~ log(Csize) * YourGroup, data = gdf, SS.type = "I", iter = 999)
anova(YourName_Full)

# Model Comparison
anova(YourName_Reduced, YourName_Full)

# Pairwise Comparisons
YourName_Pairwise <- pairwise(fit = YourName_Full, fit.null = YourName_Reduced, groups = gdf$YourGroup)
summary.pairwise(YourName_Pairwise, test.type = "dist", stat.table=TRUE)
summary.pairwise(YourName_Pairwise, test.type = "VC", stat.table=TRUE)
summary.pairwise(YourName_Pairwise, test.type = "DL", stat.table=TRUE)

# Plot Allometry with Prediction Line

# Set colors
color <- c("YourColour","YourColour","YourColour", "etc")
e<-as.factor(gdf$YourGroup)
color <- color[as.numeric(as.factor(e))]

# Plot allometry
plotAllometry(YourName_Full, size = gdf$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color, cex=1.5)

# Add legend
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))


##############################################
# HYPOTHESIS 2: SEX AND ECOLOGICAL INTERACTION #
##############################################


## Filter Data for different groups (e.g., Remove Unknown Sex and Offshore groups) ##
Indices <- filenameSlidLandmarks.dt$spec$Sex != "Unknown" & filenameSlidLandmarks.dt$spec$YourGroup != "Offshore"
filenameSlidLandmarks.dt$gpa.sh <- filenameSlidLandmarks.dt$gpa.sh[, , indices, drop = FALSE]
filenameSlidLandmarks.dt$cs <- filenameSlidLandmarks.dt$cs[indices, drop = FALSE]
filtered_spec <- filenameSlidLandmarks.dt$spec[filenameSlidLandmarks.dt$spec$Sex != "Unknown" & filenameSlidLandmarks.dt$spec$YourGroup != "Offshore", ]
filenameSlidLandmarks.dt$spec <- filtered_spec

gdfCleaned <- geomorph.data.frame(
  shape = filenameSlidLandmarks.dt$gpa.sh, 
  sex = filenameSlidLandmarks.dt$spec$Sex, 
  YourGroup = filenameSlidLandmarks.dt$spec$YourGroup, 
  ind = filenameSlidLandmarks.dt$spec$LabID, 
  Csize = filenameSlidLandmarks.dt$cs)


# Reduced Model (Common Allometry)
YourName_Reduced_Interaction <- procD.lm(shape ~ log(Csize) + YourGroup + sex, data = gdfCleaned, SS.type = "II", iter = 999)
anova(YourName_Reduced_Interaction)

# Full Model (Unique Allometry)
YourName_Full_Interaction <- procD.lm(shape ~ log(Csize) * sex * YourGroup, data = gdfCleaned, SS.type = "II", iter = 999)
anova(YourName_Full_Interaction)

# Model Comparison
anova(YourName_Reduced_Interaction, YourName_Full_Interaction)

# Pairwise Comparisons
YourName_Pairwise_Interaction <- pairwise(fit = YourName_Full_Interaction, fit.null = YourName_Reduced_Interaction, groups = interaction(gdfCleaned$YourGroup, gdfCleaned$sex))
summary.pairwise(YourName_Pairwise_Interaction, test.type = "dist", stat.table=TRUE)
summary.pairwise(YourName_Pairwise_Interaction, test.type = "VC", stat.table=TRUE)
summary.pairwise(YourName_Pairwise_Interaction, test.type = "DL", stat.table=TRUE)


# Plot Allometry
# Set color 
color2 <- c("YourColour","YourColour","YourColour", "etc")
e<-as.factor(gdfCleaned$YourGroup)
color2 <- color2[as.numeric(as.factor(e))]


Pch1 <- c(15, 17)
f <- as.factor(gdfCleaned$sex)
Pch1 <- Pch1[as.numeric(as.factor(f))]

# Plot
plotAllometry(YourName_Full_Interaction, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine", pch = Pch1, col = color2, cex=1.5)

# Add legend
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))


###################################
# HYPOTHESIS 3: SEXUAL ALLOMETRY #
###################################

# Reduced Model
YourName_Reduced_Sex <- procD.lm(shape ~ log(Csize) + sex, data = gdfCleaned, SS.type = "I", iter = 999)
anova(YourName_Reduced_Sex)

# Full Model
YourName_Full_Sex <- procD.lm(shape ~ log(Csize) * sex, data = gdfCleaned, SS.type = "I", iter = 999)
anova(Allometry.Sex)

# Model Comparison
anova(YourName_Reduced_Sex, YourName_Full_Sex)

# Plot Allometry
color3 <- c("YourColour","YourColour")
e<-as.factor(gdfCleaned$sex)
color3 <- color3[as.numeric(as.factor(e))]

plotAllometry(YourName_Full_Sex, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine", 
              pch = 16, col = color3, cex=1.5)

###############################
# COMBINE ALL PLOTS IN ONE FIGURE
###############################

par(mfrow=c(1, 3)) # 1 row, 3 columns
plotAllometry(YourName_Full, size = gdf$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color, cex=3)
plotAllometry(YourName_Full_Interaction, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color2, cex=3)
plotAllometry(YourName_Full_Sex, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color3, cex=3)


############ 3D Clouds, with Offshore ##########

# Load required libraries
library(readr)
library(ks)
library(rgl)

YourName_Full <- procD.lm(shape ~ log(Csize)*YourGroup,  
                          data = gdf, SS.type = "I", print.progress = FALSE,  iter = 999)

color <- c("YourColour","YourColour","YourColour", "etc")
e<-as.factor(gdf$YourGroup)
color <- color[as.numeric(as.factor(e))]

PlotName <- plotAllometry(YourName_Full, size = gdf$Csize, logsz = TRUE, method = "size.shape", pch = 16, col = color, cex=1.5)
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))


# Create density clouds for the first 3 PCs of the above PCA plot 
Score <- PlotName$size.shape.PCA$x[,1:3]
Labels <- factor(x=filenameSlidLandmarks.dt$spec$YourGroup)

Hscv1.GM2 <- Hkda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filenameSlidLandmarks.dt$spec$YourGroup), bw = "scv", pre = "sphere")
Tt.kda3.GM2 <- kda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filenameSlidLandmarks.dt$spec$YourGroup), Hs = Hscv1.GM2)

plot(Tt.kda3.GM2, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")

text3d(PlotName$size.shape.PCA$x[,1:3], texts = filenameSlidLandmarks.dt$spec$YourGroup, pos = 1, cex = 0.5)


############ 3D Clouds, excluding a group (e.g., Plot Coastal Only and remove the group offshore) ##########

## Filter Data (e.g, Remove Offshore Group) ##

Indice <- filenameSlidLandmarks.dt$spec$YourGroup != "Offshore" 
filenameSlidLandmarks.dt$gpa.sh <- filenameSlidLandmarks.dt$gpa.sh[, , Indice, drop = FALSE]
filenameSlidLandmarks.dt$cs <- filenameSlidLandmarks.dt$cs[Indice, drop = FALSE]
filtered_spec <- filenameSlidLandmarks.dt$spec[filenameSlidLandmarks.dt$spec$YourGroup != "Offshore", ]
filenameSlidLandmarks.dt$spec <- filtered_spec

# Updated gdf 
gdfNoOffshore <- geomorph.data.frame(shape = filenameSlidLandmarks.dt$gpa.sh, YourGroup = filenameSlidLandmarks.dt$spec$YourGroup, ind = filenameSlidLandmarks.dt$spec$LabID, Csize=filenameSlidLandmarks.dt$cs)

# Allometry
YourName_Full <- procD.lm(shape ~ log(Csize)*YourGroup,  
                                           data = gdfNoOffshore, SS.type = "I", print.progress = FALSE,  iter = 999)

# 3D plot #
color4 <- c("YourColour","YourColour","YourColour", "etc")
ee<-as.factor(gdfNoOffshore$YourGroup)
color4 <- color4[as.numeric(as.factor(ee))]

PlotName <- plotAllometry(YourName_Full, size = gdfNoOffshore$Csize, logsz = TRUE, method = "size.shape", pch = 16, col = color4, cex=1.5)
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))


# 3D CLoud plot  
Score <- PlotName$size.shape.PCA$x[,1:3]
Labels <- factor(x=dt.pca.bending$YourGroup) 

Hscv1.GM <- Hkda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filtered_spec$YourGroup), bw = "scv", pre = "sphere")
Tt.kda3.GM <- kda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filtered_spec$YourGroup), Hs = Hscv1.GM)

plot(Tt.kda3.GM, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")

text3d(PlotName$size.shape.PCA$x[,1:3], texts = filtered_spec$YourGroup, pos = 1, cex = 0.5)

A local Image A local Image